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Abstract 



Using standard calculus, explicit formulas for the one-dimensional continuous and 
discrete homotopy operators are derived. It is shown that these formulas are equiv- 
alent to those in terms of Euler operators obtained from the variational complex. 

The continuous homotopy operator automates integration by parts on the jet 
space. Its discrete analogue can be used in applications where summation by parts 
is essential. Several example illustrate the use of the homotopy operators. 



The calculus-based formulas for the homotopy operators are easy to implement in 
computer algebra systems such as Mathematica and Maple. The homotopy operators 
can be readily applied to the symbolic computation of conservation laws of nonlinear 



partial differential equations and differential-difference equations. 
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^ 1 1 Introduction 

h ; 

During the development of symbolic algorithms [4,5,6] for the computation 
of conservation laws of nonlinear partial differential equations (PDEs) and 
nonlinear differential-difference equations (DDEs) we encountered powerful 
tools from the calculus of variations and differential geometry that deserve 
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attention in their own right. This paper focuses on one of these tools: the 
homotopy operator. 

Inspired by work of Kruskal et al. [9] , we give a straightforward derivation of 
formulas for the continuous homotopy operator and its discrete counterpart. 
We show that our formulas are equivalent to those in terms of Euler operators 
obtained from the variational complex [6] . For lack of space, we only cover the 
one-dimensional cases (ID) with independent variable x or lattice variable n. 
The generalization to multiple independent variables is cumbersome [5]. 

The continuous homotopy operator goes back to Volterra's work [13] on the 
inverse problem of the calculus of variations. The homotopy operator also 
appears in the proof of the converse of Poincare's lemma [12], which states 
that exact differential forms are closed and vice versa (at least on a star-shaped 
domain in Euclidean space). Poincare's lemma is a special case of the so-called 
de Rham complex, where one investigates the equivalence of closedness and 
exactness of differential A;— forms in generality. The key to exactness proofs 
of various complexes (such as the de Rham complex) is the construction of 
suitable homotopy operators [12]. 

In basic terms, the ID continuous homotopy operator reduces the problem of 
integration by parts on the jet space to a sequence of differentiations followed 
by a single definite integration with respect to an auxiliary variable. In 2D 
and 3D, the homotopy operator allows one to invert the total divergence op- 
erator [12] . Irrespective of the number of independent or dependent variables, 
the problem can be reduced to a single definite integral. Applications of the 
continuous homotopy operator in multi-dimensions can be found in [4,5,12]. 

Likewise, the discrete homotopy operator is a tool to invert the forward dif- 
ference operator whatever the application is. It circumvents summation by 
parts by applying shifts and differentiations followed by a one-dimensional in- 
tegration with respect to an auxiliary variable [6]. Applications of the discrete 
homotopy operator are given in [5,6]. 

As shown in [8,10], the parallelism between the continuous and discrete cases 
can be made rigorous as both can be formulated in terms of variational bi- 
complexes. We do not use the abstract framework in order to make this paper 
accessible to as wide an audience as possible. Aficionados of de Rham com- 
plexes should consult [2,3] and [8,10,11]. The latter set of papers covers the 
discrete variational bicomplex. 

Several examples illustrate the inner workings of the homotopy operator at 
the calculus level, without "wedges and hooks" or differential forms. Avoid- 
ing sophisticated arguments from differential geometry, we can introduce the 
powerful concept of homotopy operators to a wider audience. 
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In [4,5,6] we apply homotopy operators to the symbolic computation of con- 
servation laws of nonlinear PDEs and DDEs. Beyond DDEs, the discrete ho- 
motopy operator is useful in the study of difference equations [8,10,11]. 

Despite their universality and applicability, homotopy operators have not been 
implemented in major computer algebra systems (CAS) like Mathematica and 
Maple. CAS offer few reliable tools for integration (or summation) by parts of 
expressions involving unknown functions and their derivatives (or shifts). We 
hope that the calculus-based formulas for the homotopy operator presented in 
this paper will lead to more sophisticated integration algorithms within CAS. 

In summary, our paper has the following objectives: (i) Give a straightforward 
derivation of the ID continuous homotopy operator, (ii) Present the discrete 
homotopy operator by analogy with the continuous case, (iii) Illustrate the 
inner workings of the homotopy operators at the calculus level, (iv) Present 
alternate, readily applicable formulas for the homotopy operators which lead 
to efficient and fast symbolic codes for integration and summation by parts. 

2 Derivation of the continuous homotopy operator in ID 

In [4,5], we presented the homotopy operator and referred to [12] for a proof, 
which involves working with differential forms. Inspired by [9], we give a 
calculus-based derivation of the homotopy operator Ti u M. For simplicity and 
clarity, we show the derivation for one dependent variable u and one indepen- 
dent variable x (henceforth referred to as the ID case). 



To do so, we first introduce a "degree" operator Ai and its inverse, the to- 
tal derivative T> x , and the Euler operator C^ x y The latter is also called the 
variational derivative or Euler-Lagrange operator. 



The calculations below are carried out in the jet space where one treats 
u,u x ,U2x, etc., as independent. As usual, u x = ^,U2 X = f^l, etc. The opera- 
tors act on f(u,u x ,u 2x , ■ ■ ■ ,um x )- Such functions are called differential func- 
tions [12]. Throughout this paper we assume that the differential functions 
lack constant terms and that the upper bounds in the summations equal the 
order M of the differential function the operators are applied to. 

Definition 1 The degree operator M. is defined by 

m ^ d f d f d f d f d f 

Mf = ^u ix - — = u — + u x - h« 2 ^ 1 \-u Mx -^ , (1) 

du ix du du x du 2x ou Mx 

where f is a differential function of order M. 

Example 2 If f = u p u q x u r Zx1 where p,q, and r are non-negative integers, then 



3 

g = Mf = Y,u lx i ± - = (p + q + r) u p u q x u r 3x . (2) 



df 



=o 



du 
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Thus, application of Ai to a monomial results in multiplication of the mono- 
mial with its degree, i.e. the total number of factors in that monomial. 

We use the homotopy concept to construct the inverse operator, M^ 1 . Given 
a differential function g(u), let g[Xu] denote g(u) where u is replaced by Aw, 
u x is replaced by Xu x , etc., where A is an auxiliary parameter. We now show 

M- 1 g(u) = Jg[M-. (3) 
o 

Indeed, if g(u) has order M then so does g[Xu], and 

d n , dg[\u) d\u ix lit dg[\u] 1 

T\ 9 M = g 1)A^-— = A g x ~dxtix~ = X ' < 4 > 



Upon integration of both sides with respect to A, we get 
/ ^[A«] rfA = g[\u]\ x x zl = g(u) - g(0) = j Mg[Xu\ ^ = mJ g[Xu] ^.(5) 



Assuming g(0) = and applying M.~ l to both sides of (5), Eq. (3) readily 
follows. The assumption g(0) = restricts the choice for /. We only con- 
sider differential functions involving monomials in u,u x , etc., and on occasion 
multiplied by sin-u or cos-u. 

Example 3 For g in (2), we have g[Xu] = (p + q + r)X p+q+r u p u q x u\ x . 
Using (3), 

l 

M- 1 g = J(p + q + r) X p+q+r ~ l u p u q x u\ x dX (6) 



= u p u q x u\ x X p+q+1 



= u p u q x ul x . (7) 



Definition 4 The total derivative operator V x is defined by 

<n f 9 f df df df 

V x f = g u {l+1)x ^- = u x - + u 2x — + ■■■ + u {M+1)X g^-. (8) 

Example 5 If f = u p u q x u 3x , then 
3 df 

V x f = u (i+-i>f, — = P uP lu t u lx + Q uPu l lu 2xU r 3x + ru p u q x u r 3x l u Ax . (9) 

i=0 OUi x 

Theorem 6 The operators M. and V x commute as do Mr x and V x . 

Proof. The proof that M. and V x commute is straightforward: applying M, 
to D x f gives the same result as applying T> x to Aif, using standard calculus 
manipulations. Proving that Ai^ 1 commutes with T> x is then immediate. ■ 
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Definition 7 The continuous Euler operator of order zero (variational deriva- 

,(0) 

-u(x) 



tive) C^L) is defined [12] by 



M QJ 



fc=0 



^-T> x — + V 2 -^--V 3 -^- + --- + (-l) M V M — (10) 
<9« ' du x x du 2x x du 3x x du Mx 



Definition 8 A differential function f of order M is called exact if there 
exists a differential function F of order M — 1 so that f = V X F. 

Theorem 9 A necessary and sufficient condition for a differential function 
f to be exact is that C^ x ^f = 0. 

Proof. A proof is given in e.g. [9]. 

Example 10 Let f—2u x u 2x cos u—u x sin u. Note that f=V x F with F— u 2 x cos u. 



We show that f is indeed exact. Using (10), we readily verify that 
mo) f _df df 2 df 



u(x)J du x du x x du 2x 

= —2u x U2x sinw — u x cosu — T> x [2u2 X cosu — 3u x sin u] +T>l[2u x cosu] 
= —2u x U2 X sin u — u x cos u — [2u 3x cos u — 8u x U2 X sin u — 3u x cos u\ 
+ [2u 3x cos u — 6u x u 2x sin u — 2u\ cos u] = . (11) 



Definition 11 The continuous homotopy operator with variable u(x) is 

n u(x) f = J (/„/) [Am] y , (12) 
o 

where the integrand I u f is given by 

M-l M a? 

U = E «<* E (-2>*) Mi+1) 7T- (13) 

Theorem 12 Given an exact differential function f of order M one has 
F = V^f = ffdx = H u(x) f. 

Proof. We multiply C^ x) f = El^o(-^) fe ^{; by u to restore the degree. 
Next, we split off ujjt. Then, we integrate by parts and split off u x J^-. We 
repeat this process until we split off umx qJ m ■ In detail, 



=«|-x>.(«g(-f.)^) + «.g(-f.)' 



<9w 



fca: 
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/ M 

V k=i 



df df ^ f ^ Qf 

. U J +U • ' ' ' n ' 



+ Ux j2(-v x ) k ~ 2 ^- + U2x j2(-v x y 

k=2 aUkx / k=2 



du, 



kx 



df df df ^ ( ™, ^ , df 



du du x du M x V t~L du kx 

^EW W #- + • • • + E (-^) fc_M | L ) 

fc=2 ""to fc=M ua kx J 

M p)f ( M-1 M p)f 

=Y,u^-v x e ^ E (-*>«) Mi+1) #- 

(M-1 M Of \ 

E^ E (-^) Mm) ^- ■ (i4) 
i=0 fc=i+l toy 

Since / is exact we have C^Kf = 0. Eq. (14) then implies 

Mf = v x [Y: u ix E (-v x ) k -^ +l) -^- . (15) 

\ i=0 k=i+l ° Ukx j 

Applying M.' 1 and using M.~ l T> x = T> X M.~ X from Theorem 6, we obtain 

(M-1 M f)f \ 

M- 1 ^^ E {-V x ) k ~ {i+1) f- ■ (16) 
i=0 fc=i+l toy 

Applying "D" 1 and using (3), we get 

/• ( M - x M df \ rl\ 

f = = j E «*, E (-^) M<+1) £r t = (17) 

g y i=o fc=«+i to y 

using (12) and (13). ■ 

Example 13 Let f = 2u x u 2x cosu — -u^sin-u which is exact as shown in Ex- 
ample 10. Using (13) with M = 2, we readily compute 

Iuf = 2_^u ix > (-£>*) ^ = uV *Vx — )+M~ — ) 

S) fc=i+l 9 ^to <9m x dw 2:c dw 2a; 

= u(2u 2x cosm — 3w^ sinw) — uV x (2u x cosm) + u x {2u x cosw) 
= — wii^ sin m + 2it^ cos m. (18) 

F — Hu(x)f — J (luf) [^ u ] -y — y ( — \ 2 uu 2 x sin(Aw) + 2Aw^ cos(Aw)) g?A 



= u 2 x cosu. (19) 
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Thus, application of (12) yields F without integration by parts with respect 
to x. Indeed, F can be computed via repeated differentiation followed by a 
one-dimensional integration with respect to an auxiliary variable A. 

3 Alternate form of the continuous homotopy operator in ID 

Formulas (12) and (13) are valid for one dependent variable u. In [5] we pre- 
sented a calculus-based formula for the homotopy operator for N dependent 
variables in ID based on work by Anderson and Olver in [12, p. 372]: 

«»(*)/ = /E(W)[Au]-, (20) 
o J =1 

where is the jth component of u = (u^\ u^ 2 \ • • • , • • • , u^). The 
integrand, 

mW-i 

W = E 0j(« w 4$o/)> ( 21 ) 

i=0 

where is the order of the variable in /, involves the continuous ID 
higher Euler operators [9,12] defined as follows. 

Definition 14 The continuous higher Euler operators for component u^\x) 
are 

^ w / = f 0<-*.>"aJb (22 > 

where (?) is the binomial coefficient. 

Note that the higher Euler operator for i = and one dependent variable 
= u(x) matches the variational derivative (10). 

In the case of one dependent variable, u, we denote by M. After substi- 
tution of into (21), we obtain 

"-t*('Mi+i) i -*^&)- (23) 

Theorem 15 The integrands (13) and (23) are equal. 

Proof. Starting from (23), we use Leibniz's rule to propagate V x to the right 

df 



M-l 




(24) 
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Next, we interchange the sums over % and j (to bring Uj X up front), followed 
by an interchange of the sums over i and k (to bring V x and df/duk x outside 
the sum over i). So, 

M-1 M pif 

= E u 3X e (-^r u+1) ir^ (25) 

j=0 k=j+l UUkx 

where we have used the identity 

which is straightforward to prove using mathematical induction. ■ 

Homotopy operator (20) can easily be implemented in CAS. In our experience, 
integrand (13) leads to a more efficient and faster algorithm for it requires 
substantially less differentiations then (23). 

4 Discrete Euler and homotopy operators 

We now turn to the discrete analogues of differential functions, Euler operators 
and homotopy operators. For simplicity we consider the case of one lattice 
variable n which results, for example, from a discretization of the variable x. 
We allow iV dependent variables, i.e. u n = (u£\ vf£\ • • • , ■ ■ ■ , u^). For 
simplicity, in the examples we denote these components by u n ,v n , etc. 

By analogy with V x and V' 1 , we define shift operators acting on f n (u n ). 

Definition 16 D is the upshift operator (also known as forward- or right- 
shift) such that D f n = f n +i- Its inverse, D -1 , is the down-shift operator (or 
backward- or left-shift) such that D _1 f n = / n _i. The identity operator is 
denoted by I. Lastly, A = D — I is the forward difference operator so that 
A / n = (D- I) /„ = /„+!-/„. 

Given are functions /„ in discrete variables u n ,v n , . . . and their up and down 
shifts. If f n (u n - p , v n _ r , ■ ■ ■ , u n , v n , ■ ■ ■ , u n+q , v n+s ) with p ^> r involves negative 
shifts, one must first remove these by replacing f n by f n = D p / n . From this 
point we assume that all negative shifts have been removed. 

Definition 17 /„ is called exact if it is a total difference, i.e. there exists a 
F n so that f n = AF n . 
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Example 18 Let 

f n = sin(w n+3 ) cos 2 (w 2 +2 ) - sin(w n+ i) cos 2 (v 2 ). (27) 
By hand, we readily verify that f n = AF n with 

F n = sin(w n+2 ) cos 2 (w 2 +1 ) + sin(w n+1 ) cos 2 (w 2 ). (28) 
So, f n is exact. 

Below we address the following questions, (i) How can one test whether or not 
/„ is exact? Equivalently, how does one know that F n exists in closed form? 
(ii) Can one compute F n = A~ 1 f n in a way analogous to the continuous case? 

Definition 19 £ ( % is the discrete Euler operator of order zero ( discrete vari- 
ational derivative) for component defined [1] by 



U "' , , <).(„ _d 

k=0 ° U n+k ° Un k=0 



C { %f n = V D- k ^- = 4,V D~ k f n 



= (I + D- 1 + D- 2 + ■ • ■ + D-" W) ) /„, (29) 

OUn 

where is the highest shift of u$ occurring in f n . 

With respect to the existence of F n , the following exactness criterion is well- 
known and frequently used [1]. 

Theorem 20 A necessary and sufficient condition for a function f n with pos- 
itive shifts to be exact is that C^ f) f n = 0, j — 1, 2, 3, • • • , N. 

Proof. A proof is given in e.g. [7]. ■ 

Example 21 To test that (27) is exact we apply (29) to f n for each component 
of Un = (iiW,4 2 ') = (u n ,v n ) separately. For component u n with maximum 
shift M« =3 we readily verify that ££)/„ = ^ (I + D" 1 + D~ 2 + D~ 3 ) f n = 

0. Similarly, for component v n with maximum shift = 2 we check that 
4°J/n = £(I + D- 1 + D- 2 )/„ = 0. 

Next, we compute F n so that f n = AF n = F n+1 — F n . 
Definition 22 The discrete homotopy operator for u„ is 

In j\ 
-HuJ n = /Efc/jMy (30) 
o i =1 

where the integrand I a) f n is given by 

M<')-1 o M<i) 

I&f*= E V&Af E D~ (fc " l) / n . (31) 
i=0 OU n+i k =i+l 
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As in the continuous case, [I <j)f n ) [AuJ means that after / y)/ n is computed 

one replaces u n by Au n , u n+1 by Au n+1 , etc. 

One can use the following theorem [6,8,10] to compute F n . 

Theorem 23 Given an exact function f n one has F n = A _1 /n — 1~tu n fn- 

Proof. The proof is similar to that of Theorem 12. For simplicity we show 
it for one dependent variable, u$ = u n , and we denote by M. After 

multiplication of C^f n = Sjtio D -fe /„ with u n , we isolate u n ^- and apply 
A to the remaining term in order to sum by parts. Next, we isolate u n+ i gj n +1 , 
followed by another summation by parts. 
This process is repeated M — 1 times, so that 

M /m-i a m \ 

UnJC^f n = E u n+i ^f- - A £ u n+l -^- £ V-^f n . (32) 
i=o ua n+t y j =0 <Ja n+t k=i+1 J 

The inverse of the discrete operator M.f n = Y^Lq u n+i g~ is computed as 
M.~ x f n = Jo 1 fn[^u n ] Since / n is exact and 1 and A commute, we get 

' /M - 1 d " \ r , ,dA 



F n = A- 1 /™ = / E Vio E D" ( " _i) /n [Au„] -T- = n u J n (33) 

using (30) and (31). 



Thus, the homotopy operator reduces summation by parts needed for the 
inversion of A to a set of shifts and differentiations followed by a single definite 
integral with respect to a scaling parameter A. 

Example 24 We return to (27) where f n involves u n and v n with maximal 
shifts M« = 3 and = 2. Using (31), we get 

I u Jn = u n —^D- 1 + D- 2 + D- 3 )f n +u n+1 - (D- 1 +D- 2 )/ n +w n+2 - D-7„ 

OU n OU n+ i OU n+2 

= Un du~ ( sin ( Mn + 2 ) cos2 ( v n+i) + sin(w n+ i) cos 2 (v 2 ) 
- sin(w n _i) cos 2 (t> 2 _ 2 ) - sm(u n - 2 ) cos 2 (v 2 „ 3 
+«n+i o d (sin(w n+2 ) cos 2 (w 2 +1 ) - sm(u n ) cos 2 (w 2 _ 1 ) 
+ sin(M„ +1 ) cos 2 (v 2 ) - sin(w„_i) cos 2 (w 2 _ 2 )) 
+u n+2 ^ (sin(w n+2 ) cos 2 (w 2 +1 ) - sin(u n ) cos 2 (w 2 _ 1 )) . 

OU n+ 2 x ' 

= u n+ i cos(w n+ i) cos 2 (w 2 ) + u n+2 cos(w n+2 ) cos 2 (w 2 +1 ). (34) 
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and, analogously, 

J«Jn = -4 (yl sin(u n+ i) cos(^) sin(^)+^ +1 sin(M„ +2 ) cos(v 2 +1 ) sin(^ +1 )).(35) 
Based on (30), we compute 



{Iujn + Iv n fn) [AU„] y 




1 



= J («n+l COs(AM n+ i) COS 2 (A 2 ^) + W n+2 COs(A«„+ 2 ) COS 2 (A 2 V 2 +1 ) 





— 4 A v 2 sin(Au n+ i) cos(A 2 t> 2 ) sin(A 2 t> 2 
-4Aw 2 +1 sin(AM„ +2 )cos(A 2 t; 2 +1 )sin(A 2 t; 2 +1 )) d\ 
= sin(w n+2 ) cos 2 (v 2 +1 ) + sm(u n+1 ) cos 2 (v 2 ), (36) 

which agrees with (28) previously computed by hand in Example 18. 

5 Alternate form of the discrete homotopy operator 

In [5] we presented the following formula for the discrete homotopy operator 
with one lattice variable (n) and with iV dependent variables uff: 

(l u U)fn) [Au n ] ^ 

i =1 



n u Jn = /E(^)/n)[Au n ] T (37) 



with 



M"0)-l 

'«>/» = E A j kW/n , (38) 



i=0 " ' 

where the discrete higher Euler operators are defined as follows. 
Definition 25 The discrete higher Euler operators for component u$ are 

<S>/« = E (*V;rro- = A*f ( -V"'/.. (39) 



Theorem 26 T/ie integrands (31) and (38) are equal, i.e. 
M«)-l / o m«) / 7 \ 



MW-l . g MO) 
=0 CW n+i fc=i+i 



= E"4i"E D- (fe -^ (40) 



Proof. The proof is analogous to that of Theorem 15. ■ 
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6 Conclusions 

In this paper we derived formulas for the ID continuous and discrete homotopy 
operators. We showed that our calculus-based formulas are equivalent to those 
obtained from the variational complexes. The simplified formulas no longer 
involve higher Euler operators which makes them easier to implement and 
faster to execute in major computer algebra systems. Simplified versions of 
the continuous homotopy operator in 2D and 3D will be presented elsewhere. 
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